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Abstract 

QCD thermodynamics is considered using Wilson fermions in the fixed scale approach. The 
temperature dependence of the renormalized chiral condensate, quark number susceptibility 
and Polyakov loop is measured at four lattice spacings allowing for a controlled continuum 
limit. The light quark masses are fixed to heavier than physical values in this first study. 
Finite volume effects are ensured to be negligible by using approriately large box sizes. The 
final continuum results are compared with staggered fermion simulations performed in the fixed 
Nt approach. The same continuum renormalization conditions are used in both approaches and 
the final results agree perfectly. 
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1 Introduction 



As the early universe evolved a transition occured at temperatures T sa 150 — 200 MeV, which is 
related to the spontaneous breaking of chiral symmetry in QCD. The nature of the QCD transition 
[I] affects our understanding of the history of the universe; see e.g. [2J. 

Extensive experimental work is currently being done with heavy ion collisions to study the QCD 
transition, most recently at the Relativistic Heavy Ion Collider, RHIC and at the Large Hadron 
Collider, LHC. Both for the cosmological transition and for RHIC/LHC, the net baryon densities 
are quite small, thus the baryonic chemical potentials fi are much less than the typical hadron 
masses, /j, is below 50 MeV at RHIC, even smaller at LHC and negligible in the early universe. 
Thus, a calculation at fj, = is directly applicable to the cosmological transition and most probably 
also determines the nature and absolute temperature of the transition at RHIC/LHC. Therefore, 
we carry out our analysis at /i = 0. Given the far-reaching implications, it is desirable to perform 
this calculations in a framework that is conceptually clean. For a pedagogical review on the QCD 
transition at /i = and \i > see e.g. [3] and for the first continuum result at /x > - namely, for 
the curvature on the chemical potential vs. temperature plane, see [3]. 

When we analyze the absolute scale or any other question related to the T > QCD transition 
for the physically relevant case two ingredients are quite important. 

First of all, one should use physical quark masses. The nature of the transition strongly depends 
on the quark mass. Lattice studies and effective models showed that in the three flavor theory for 
small or large quark masses the transition is a first order phase transition, whereas for intermediate 
quark masses it is an analytic crossover. Since the nature of the transition influences the absolute 
scale T c of the transition - its value, mass dependence, uniqueness etc. - the use of physical quark 
masses is essential for the determination of T c , too. The absolute scale then goes into all observables. 
Whereas it is relatively easy to reach the physical value of the strange quark mass m s in present 
day lattice simulations, it is much more difficult to work with physical up and down quark masses 
TTt u( i , because they are much smaller: m a /m u d ~ 28. In calculations with m s /m u d smaller than 28 
the strange quark mass is usually tuned to its approximate physical value, whereas the average up 
and down quark masses are larger than the physical value. 

Secondly, the nature and other characteristics of the T > QCD transition are known to suffer 
from discretization errors [SJ|B]. Let us mention one example which underlines the importance of 
removing these discretization effects by performing a controlled continuum extrapolation. The three 
flavor theory with a large, a 0.3 fm lattice spacing and standard staggered action predicts a critical 
pseudoscalar mass of about 300 MeV [7]. This point separates the first order and cross-over regions. 
If we took another discretization, with another discretization error, the critical pseudoscalar mass 
turns out to be much smaller, well below the physical pion mass of 135 MeV. The only way to 
determine the physical features of the transition is to carry out a careful continuum limit analysis. 
It can be safely done only in the so-called scaling regime. This regime is reached when the lattice 
spacing a is sufficiently small, smaller than some a max . Dimensionless combinations of observables 
approach their continuum limit value (within their error bars) in the scaling regime with a correction 
term c ■ a n . Here c, n and a max depend on the action and on the dimensionless combination. The 
values of c and a ma x are typically unknown, whereas the form of the action and the observables 
provide the value for n, usually without performing any simulations. To carry out a controlled 
continuum extrapolation at least three lattice spacings in the scaling regime are needed. Two 
points will always lie on a two parameter c • a n curve, independently whether the lattice spacings 
are smaller than a max or not; the third point indicates if one reached the scaling regime. 

It is numerically very demanding to fulfill both conditions. There are only a few cases, for which 
this has been achieved. Within the staggered formalism there are full results such as the nature of 
the transition [I] , the transition temperature [8HT0] , equation of state [9] and fluctuations [11] . 

At fi > lattice computations are more costly, see e.g. [12]. However the curvature of the 
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phase line separating the hadronic and quark-gluon phases is already determined in the continuum 
limit [3]. Results for the possible critical point only simulations on coarse lattices are available; see 

e.g. na. 

It is important to note that fulfilling the second condition without fulfilling the first one still 
leads to universal results. In other words continuum extrapolated results with non-physical quark 
masses are universal. Independently of the action, simulation algorithm, scale setting procedure, 
they provide the same answer once the quark mass is fixed which is a non-trivial issue, but can 
be done e.g. by fixing the pion to Omega and kaon to Omega mass ratios: M^/Mq and Mk/Mq. 
These results are not the same as they are for physical quark masses, but they are well defined and 
unique. Contrary to this universality, fulfilling the first condition (physical quark mass) but not the 
second one (continuum extrapolation) leads to non-universal, non-physical results. These results 
still have unknown discretization errors. 

Once the available computational resources are not enough to fulfill both conditions it is more 
advisable to carry out calculations with non-physical quark masses but perform the continuum 
limit extrapolation. As we have seen such results are universal and can be cross-checked with other 
results obtained by other fermion formalisms, actions etc. For some recent Wilson thermodynamics 
results see (I3H20]. 

In this paper we determine the temperature dependencies of a couple of obscrvablcs (chiral 
condensate, strange susceptibility, Polyakov loop) in 2+1 flavor QCD. We use Wilson fermions with 
six steps of stout smearing and tree level clover improvement in the quark sector and a tree level 
improved action in the gauge sector; for the details of the action see [21]. Our pion is non-physical, 
its mass is about 545 MeV, see later for a detailed discussion for the mass. 

The structure of the paper can be summarized as follows. After this brief introductory Section 
Q] a discussion on the advantages and disadvantages of Wilson thermodynamics is presented in 
Section [5] The main features of the action and run parameters are listed in Section [31 Our choice of 
renormalization procedures for the various measured quantities are summarized in Section [4[ The 
results are given in Section [5] In Section [6] we summarize and provide an outlook. 

2 Choice of the fermion formalism 

This paper presents the first of a series, which deals with lattice QCD thermodynamics using Wilson 
fermions. The final goal is to describe several bulk observables as a function of the temperature all 
the way to the continuum limit with physical quark masses or in other words by approaching pion 
masses of 135 MeV and kaon masses of 495 MeV. While the framework has the advantage of the 
sound conceptual status of the Wilson fermion formulation (a), and we are already able to reach the 
continuum limit (6), the present work suffers from the large computational cost (c) and is therefore 
confined to larger than physical pion masses (cf). Let us discuss these aspects in some detail. 

ad a. The vast majority of the large scale lattice QCD thermodynamics projects has been carried 
out with staggered fermions. Staggered fermions are the cheapest formulation of lattice QCD and 
working with them turned out to be quite successful and provided many interesting results both at 
T = and T > 0. In their original form they describe four degenerate fermions - tastes - and one 
has to take the square root or the fourth root of the fermion determinant to describe two fermions or 
one fermion, respectively. This procedure is somewhat unattractive and there has been an ongoing 
debate in the literature whether it leads to the proper universality class. Wilson fermions do not 
raise such theoretical questions. Furthermore, staggered fermions suffer from the so-called taste 
symmetry violation. It means that instead of the physical pseudo-Goldston bosons there is a tower 
of pscudoscalars whose masses are typically well beyond the physical pion or kaon masses. The 
physical spectrum is expected to be restored only in the continuum limit. In the last five years 
some of the authors of the present paper have carried out a large scale staggered thermodynamics 
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program, nevertheless it is desirable to cross check those results by repeating the analyses with 
Wilson fermions. Its theoretical cleanness is the most important reason to perform a systematic 
study with Wilson fermions. 

ad b. As we discussed earlier lattice results are unambiguous only in the continuum limit. To 
this end we carry out our analysis at four lattice spacings. As we will see the results scale quite 
nicely and the continuum behavior can be extracted. The results are in good agreement with the 
results obtained with another fermion discretization, stout smeared staggered fermions. There is, 
however, one important conceptual difference between Wilson and staggered thermodynamics. In 
the staggered case the standard procedure is to take a given temporal extent N t to control the 
lattice spacing and use various gauge couplings, light and strange quark masses /3, m u d and m s to 
change the temperature. The set of parameters /?, m u d and m s , for which the physical content is the 
same (e.g. the ratio of the kaon decay constant and the pseudoscalar masses: fx/M^ and Jk/Mk) 
are called lines of constant physics (LCP). The fixed scale approach of Wilson thermodynamics [T3] 
is in some sense the opposite. One takes a given (3 to fix the lattice spacing and uses several Nt 
values to scan the temperature. The reason for this choice is the difficulty with the additive mass 
renormalization within the Wilson formalism. This additive term makes it particularly difficult to 
give the LCPs with good numerical precision. The fixed scale approach ensures by definition that 
the physical content remains the same (same bare parameters) even if we change the temperature. 

ad c. Wilson fermions arc usually more expensive than staggered fermions. There are at least two 
reasons for that. First of all the basic calculational step, fermion matrix x spinor multiplication, 
has four times as many floating point operations for Wilson fermions as for staggered fermions. 
Secondly, the computational costs for the inversion of the fermion matrix strongly depends on its 
condition number. The condition number for the fermion matrix has a strict bound given by the 
quark mass for the staggered case, whereas for the same quark mass with Wilson fermions the 
condition number can be much larger. It depends not only on the mass but also on the gauge 
configuration and thus fluctuates more. 

ad d. Today, the large computational costs allow to study only systems with larger than physical 
pion masses; larger mass means smaller condition number, thus smaller computational costs. As 
we discussed earlier it is more reasonable to carry out a study with larger than physical mass 
and extrapolate to the continuum limit than use smaller or even physical quark masses at one or 
two lattice spacings only. The reason for that is simple: continuum extrapolated results can be 
compared with those of other groups and/or lattice actions. A non-continuum result has still an 
unknown uncertainty due to cutoff effects. As we mentioned we compare our findings in this paper 
using Wilson fermions with our earlier results using staggered fermions. A nice agreement is found. 



3 Simulation points and techniques 

In this section the details of the Wilson simulations are outlined. We also performed staggered 
simulations in order to compare the continuum results of the two formulations. The techniques of 
the staggered simulations follow [5112] and are shortly summarized too. 



3.1 Action 

The gauge action used for the calculations was the Symanzik tree level improved gauge action [^1125] 
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Table 1: Simulation parameters. The Nt values used for the finite temperature runs and the values 
used for the zero temperature runs are separated by a comma. 

with the parameters Co = 5/3 and c\ = —1/12. The action for the femionic sector was the clover 
improved [21] Wilson action 

= az. y: E *. > ( 2 ) 

where Sp is the Wilson fermion action. Six steps of stout smearing [25] with smearing parameter 
q = 0.11 were used. The clover coefficient was set to its tree level value, csw = 1-0, which, for this 
type of smeared fermions, essentially leads to an 0(a) improved action |26j with improved chiral 
properties [2"7] . 

3.2 Simulation algorithm 

The bare masses of the u and d quarks were taken to be degenerate, therefore the configurations 
were generated using an Nf = 2 + 1 flavor algorithm. The light quarks were implemented via the 
Hybrid Monte Carlo (HMC) algorithm [35], whereas the strange quark was implemented using the 
Rational Hybrid Monte Carlo (RHMC) algorithm [223 - I n order to speed up the molecular dynamics 
calculations, the Sexton- Weingarten multiple time-scale integration scheme [30] combined with the 
Omelyan integrator [21] was employed. When all four extents of the lattice were even, the usage 
of even-odd preconditioning [32] gave an additional speed up factor of 2. For further details on the 
algorithm see |21| . 

3.3 Simulation points 

The calculations were performed at four different gauge couplings, (3 = 3.30, 3.57, 3.70 and 3.85. 
Only dimensionless ratios are measured and every dimensionful quantity is made dimensionless 
by appropriate powers of mo or m ff . The lattice spacing is also set by the physical value m,Q = 
1672 MeV and the four gauge couplings correspond to a = 0.139(1), 0.093(1), 0.070(1) and 0.057(1) 
fm, respectively. At all four gauge couplings the bare masses are tuned such that m^/mQ = 0.326(4) 
and mif/mn = 0.366(4) are constant, which means that « 545 MeV and mx ~ 614 MeV. Wc 
tried to tune the m s /m u d ratio to the same value m s /m u d = 1.5 as used in the staggered reference 
runs. We achieved this within 2% accuracy: the ratio (2rn K — m^)/m^ which in leading order chiral 
perturbation theory equals m s /m uc i is 1.530(7) for all four gauge couplings. All four values of m s 
are such that if m uc i is lowered to the physical point mx also becomes the physical kaon mass. The 
bare quark masses, spatial and temporal lattice extents are shown in table [T] while the measured 
masses are shown in tabled] In all four cases m T L s > 8. At each finite temperature point around 
1000 equilibrated trajectories were generated while around 500 at zero temperature points. 
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a [fm] 


Z A 


3.30 


0.332(3) 


0.373(3) 


1.529(2) 


0.0428(2) 


1.16(1) 


0.139(1) 


0.892(7) 


3.57 


0.319(6) 


0.359(4) 


1.531(2) 


0.02649(4) 


0.777(9) 


0.093(1) 


0.951(2) 


3.70 


0.326(5) 


0.369(5) 


1.531(3) 


0.01994(4) 


0.586(8) 


0.070(1) 


0.966(2) 


3.85 


0.314(7) 


0.358(6) 


1.528(4) 


0.01559(2) 


0.480(8) 


0.057(1) 


0.976(5) 



Table 2: Spectroscopy and Za renormalization constant results from zero temperature simulations. 
The m s /m u d column refers to (2m|- — m 2 )/m 2 . The lattice spacings are set by mo = 1672 MeV. 

3.4 Staggered simulations 

The goal of the staggered simulations detailed below was to provide a basis of comparison for the 
Wilson data. We used the staggered action with two steps of stout smearings (g = 0.15) as in 
most our thermodynamics studies e.g. [5]. The line of constant physics used there was extended 
to the fine lattices in [9 , 33 , 34 . We used the bare strange masses for each given beta as were 
determined there. The light quark mass was simply set to 2/3 of the strange mass to achieve the 
desired pion/kaon mass ratio. The lattice spacing in these staggered thermodynamics papers was 
set through the kaon decay constant. For each used gauge coupling we determined the scale directly 
with zero temperature runs at the light/strange mass ratio 2/3. Using the most precise scale setting 
observable we had access to, the wo scale [35], we checked that our original scale function was still 
correct. We proceeded to compute the mass ratios m n /m^i ~ 0.32 and Mjf/mfi — 0.36 that turned 
out to agree with the Wilson values at the level of 3% in both cases. 

For the renormalization we used dedicated zero temperature runs with parameters precisely 
matching those in the finite temperature simulations. We determined the vacuum condensate (^V)o 
and the static potential V(r), which we needed for obtaining the multiplicative renormalization of 
the Polyakov loop: Z = cxp(V(x)/2). The obtained finite renormalized Polyakov loop is then 
further transformed to our scheme (see section I4.4[) by a finite renormalization factor. One can 
select any physical scale x to remove all divergences. Since wq is our most accurately known scale, 
we used x = \/8wo, which is approximately 0.5 fm. We found that in the continuum, at this mass 
ratio we have m^WQ = 1.466(15). 

Our finite temperature staggered simulations were performed on 18 3 x 6, 24 3 x 8, 32 3 x 10 and 
36 3 x 12 lattices. After performing the necessary steps of renormalization a continuum extrapolation 
was carried out. For this extrapolation we used a cubic spline interpolation (using roughly every 
second temperatures as node points). Then for every temperature we evaluated several possible ex- 
trapolations (linear in a 2 ). The choices included using all four lattice spacings or just three of them, 
and also making the extrapolation for the reciprocal observable. The width of the weighted his- 
togram of the continuum extrapolations define a systematic error, which we added to the statistical 
error in quadrature. The final continuum extrapolation includes an additional overall percent-level 
error in the temperature axis. 

4 Renormalization 

The bare chiral condensate is divergent in the continuum limit and both additive (power-like) and 
multiplicative (logarithmic) divergences need to be removed. The resulting finite chiral conden- 
sate in the continuum limit can be compared with results obtained using other regularizations 
for instance the staggered formulation since finite continuum quantities should not depend on the 
regulator. 

The additive and multiplicative renormalizations are treated separately in the following two 
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subsections. We follow [35] which is based on [37]. The measurement of the finite renormalization 
constant Za is outlined in subsection 14.31 The Polyakov loop also needs to be renormalized and 
our scheme is defined below in subsection 14.41 



4.1 Additive renormalization of the condensate 

On dimensional grounds the bare chiral condensate contains additive divergences of the type 

(VWo) 



N 



f 



c + ci(m - m c ) + c 2 {m - m c ) +... (3) 



where mo is the dimensionful bare mass, m c is its critical value and Co is cubically, c\ is quadratically 
and C2 is linearly divergent. The coefficients Ci do not depend on the temperature hence they cancel 
in the quantity 

A^(r) = (^o)(r)-(^ )(r = o) 

N f N f U 

just like in the staggered case. The cancellation is of order 0(a~ 3 ). 

However the 0(a -3 ) term in (UJ can be explicitly removed and the following quantity is free of 
cubic divergences [SB] . 

- c = 2m P cAcZ A J d 4 x(P (x)P (0)) + ■ • • , (5) 

but of course quadratically and linearly divergent pieces are still present. Here Pq{x) is the bare 
pseudo-scalar condensate, mpcAC is the PCAC mass and Z A (go) is a finite renormalization constant 
[36]. Our conventions for the definition of Po(x) and (ipoipo) are the same as in [36]. Hence for the 
subtracted condensate we have, 



N f 

where the short hand notation 



2m PC AcZA^pp(T) + • • • , (6) 



App(T) = / d 4 x(P (x)P (0)}(T) - / d 4 x(P (x)P (0))(T = 0) (7) 



was introduced. In (J7J the cancellation between the finite and zero temperature terms is only 
0(a~ 2 ) however all additive divergences are still removed. 



4.2 Multiplicative renormalization of the condensate 

Multiplicative logarithmic divergences arc still present in both (Q} and © . Multiplying both expres- 
sions by first Zp and then by a renormalized mass mp = ttipcacZa/Zp leads to a renormalization 
group invariant quantity mp(iptjj}p. Hence using ^ and ([7J we arrive at the expressions, 

m R {^) R (T) - 2N f m 2 p CAC Z A A P p(T) (8) 
mn(H->)R{T) = m PC AcZA^4,(T) + ••■ , (9) 

where the former may be taken as the definition of (ipip)p(T) and the latter agrees with it in the 
continuum limit. Both are of course also finite. The finite cut-off corrections for ([5]) is 0(a 2 ) at tree 
level provided the action is O(a)-improvod, because the tree level improvement factors for mp and 
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{^ip)R cancel in the product. The full non-perturbative improvement of mp and (V'V)ii we expect to 
be close to tree level improvement because our smeared action with tree level improvement coefficient 
csw = 1-0 is very close to being non-perturbatively 0(a)-improved [251127] . The quantities A^^(T) 
and App(T) on the right hand sides are easy to measure on the lattice and knowing mpcAC and 
Za from zero temperature simulations allows one to define the finite and renormalization group 
invariant continuum quantity mp(rpip) p{T) in two different ways. In the continuum limit the two 
definitions should agree within errors. 

On the other hand, using the first expression in (0) to solve for tupcacZa and substituting it 
into the second expression leads to, 

A 2 7 , (T) 

m R lM) R (T) = 1 , ; + • • • , (10) 

which can directly be measured from the bare quantities A^(T) and App(T) without any need for 
measuring mpcAC or the Za renormalization factor. In [55] the above ratio was used to define the 
chiral condensate, in the present study we will use (jSJ) because we have found that this definition 
of the renormalized condensate scales best among the three possible choices (J5J), (JSJ) and (TTUl) . 

4.3 Computation of Za 

We compute the renormalization constant of the local axial vector current Za along the lines 
presented in [5211301 ■ We generated independent sets of Nf = 3 ensembles at four quark masses in 
the range m s /3 < m q < m s and approximately equal physical volume V ~ (2 fm) 4 at every /3. 
In a first step, we compute the local vector current renormalization constant Zy from the ratio 

(P (N t /2)V (t)P (0)) 
?W (P (JVt/2)Po(0)> 1 ' 

where 

P (t) = ( d 3 xP (t,x) , V Q (t) =[ d?x (^170^1 )(t,x) , (12) 



are the zero 3-momentum projected bare pseudo-scalar and vector densities. With tree level im- 
provement one has [¥TlEf2"l 

Z V (J3, m)(l + am) = (£(h) - {(h))- 1 for < h < N t /2 < t 2 < N t (13) 

We obtain Za by using the standard RI-MOM procedure [55] with the improvement technique 
of [33] to determine the ratio 

Za{P,™) = Ty{p) 

Z v (p,m) T A (p) [ ' 

from the off-shell amputated Greens functions Tr(p). The dependence of the ratio (TT51) on the 
external quark momenta p is very mild and enters into our estimate of the systematic error. We 
linearly extrapolate the resulting Za{P, m) to m = to obtain the renormalization constant Za{/3), 
see table [2] 

4.4 Polyakov loop 

The real part of the bare Polyakov loop also needs to be renormalized [5J in order to have a quantity 
with a finite continuum limit. Since there is an additive divergence in the free energy, a convenient 
choice of renormalization prescription is demanding a fixed value for the renormalized Polyakov 
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staggered 




Figure 1: Renormalized chiral condensate. The top left panel shows the Wilson results using 
the definition © at the four lattice spacings corresponding to four (3 values. The continuum 
extrapolated result is also shown by the solid band. The top right panel shows the staggered results 
in the fixed- N t approach also together with the continuum extrapolated result. The bottom panel 
compares the two continuum results. 
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loop at a fixed but arbitrary temperature T» > T c . Then the renormalized Polyakov loop Lr is 
given by 



in terms of the bare Polyakov loop Lq{T). We choose T* = 0.143ma and L* = 1.2. Other choices 
would simply correspond to other renormalization schemes. 

We imposed the same renormalization condition on the central value of our continuum extrap- 
olated staggered data so they can be meaningfully compared. 

5 Results 

We have measured three quantities at each lattice spacing. The renormalized chiral condensate is 
sensitive to the remnant of the chiral transition whereas the renormalized Polyakov loop and the 
strange quark number susceptibility are sensitive to the remnant of the confinement-deconfinement 
transition. 

Each quantity is renormalized properly so that in the continuum limit finite and regularization 
scheme independent values are obtained. The Wilson continuum extrapolation is based on a cubic 
spline interpolation to temperatures not reachable by the discrete range of N t and extrapolation 
a — > 0. We perform global fits to all our data points including lattice spacing dependence in the 
fit parameters as follows. Spline node points are randomly distributed along the temperature axis. 
The parameters of our fit are the values of the observable at the node points together with their 
lattice spacing dependence, i.e. on + c(a)c>2i for the i th node point. For the leading correction c(a) 
we used two choices, c(a) = a 2 and c(a) = aa. We had 4-6 nodepoints with 1000 different random 
node point sets each, resulting in 2 • 3 • 1000 = 6000 fits alltogether. The results of the fits are 
weighted with the fit qualities. For each temperature the median of these fits is used as our final 
result and the systematic error comes from the central 68% of the distribution and the statistical 
error from a jackknife analysis. The simulations were performed at 4 lattice spacings but not every 
observable is measured in the full temperature range. For the temperature range where our final 
results are presented we had at least 3 lattice spacings available in other words at least 3 lattice 
spacings are used for the continuum extrapolation. The obtained continuum results can then be 
compared with results from other regularizations such as the staggered formulation. 

Figure[T]shows the renormalized chiral condensate. The left panel contains the Wilson simulation 
results at four lattice spacings together with the continuum extrapolated result. The right panel 
shows the chiral condensate result for the staggered calculation at 4 fixed N t values 6, 8, 10 and 
12 corresponding to four lattice spacings. Cubic spline is used to interpolate at fixed N t and these 
are then extrapolated to the continuum. The middle panel shows the two continuum extrapolated 
results. They are in perfect agreement. 

In an earlier exploratory work we used the definition (|10[) for the renormalized chiral condensate 
resulting in much larger cut-off effects [55] • The reason for using J5J this time is exactly because 
the cut-off effects are smaller as is visible from figure Q] compared to figure 1 in [55] . 

The strange quark number susceptibility \s = T/V d 2 log(Z)/dfi 2 is a sum of two contributions, 
the connected and disconnected terms. The disconnected part is a very noisy quantity (as usual) 
and a large number of random vectors, 1200, were needed in order to evaluate it precisely. It is 
advantageous to generate the random vectors in the 12 diagonal spin-color blocks separately |15j . 
Cut-off effects may be reduced by tree level improvement. The measured susceptibility is divided 
by its infinite volume, massless Stefan-Boltzmann limit at the given finite Nt temporal extent. The 
continuum limit is unchanged by this improvement since the Stefan-Boltzmann limits are 1 for 
Nt — s- oo; see table [5] 




(15) 
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Wilson 



staggered 





Figure 2: Strange quark number susceptibility. The top left panel shows the Wilson results at 
the four lattice spacings corresponding to four (3 values. The continuum extrapolated result is 
also shown by the solid band. The top right panel shows the staggered results in the fixed- N t 
approach also together with the continuum extrapolated result. The bottom panel compares the 
two continuum results. 
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N t 


8 


10 


12 


14 


16 


18 


20 


22 


24 


26 


SB 


1.522 


1.265 


1.161 


1.110 


1.081 


1.062 


1.049 


1.040 


1.034 


1.028 



Table 3: The tree level improvement factors of the strange quark number susceptibility for Wilson 
fcrmions. The values shown are the free, infinite volume, massless Stefan-Boltzmann limits at given 
Nt- For ft = 3.30 and 3.57 only N t > 8, for P = 3.70 only Nt > 12 and for P = 3.85 only N t > 14 
is used. 

The results for the four lattice spacings are shown on figure [5] together with the staggered results. 
What is plotted in both cases is the aforementioned tree level improved quark number susceptibility 

The number of data points is less for the quark number susceptibility than for the chiral con- 
densate because for the condensate odd Nt values are used as well. In such a setup even-odd 
preconditioning can not be used resulting in a much slower simulation and measurement. The large 
number of random vectors needed also forced us to also skip some even Nt values in some cases 
only having those that are divisible by 4. 

Discretization errors are comparable to the chiral condensate. The results are shown on figure [5] 
again with the Wilson results on the left, staggered results on the right and the comparison of the 
two continuum results in the middle. Again the agreement between the two continuum extrapolated 
results is perfect. 

The third quantity we measured is the rcnormalizcd Polyakov loop which is sensitive to the 
confincmcnt-deconfincmcnt transition similarly to the quark number susceptibility. The results 
are shown on figure [3] in the same format as before; Wilson results on the left, staggered results 
on the right with the comparison of the two continuum extrapolated curves in the middle. The 
bare Polyakov loop can be measured quite precisely but renormalization in our scheme at low 
temperatures causes the errors to increase because of the division and high power in equation (1151) . 
For higher temperatures the errors are smaller where the renormalization condition (|15p does not 
introduce large uncertainties. 

6 Summary and outlook 

In this work finite temperature QCD was studied using Wilson discretization for 2 + 1 flavors of 
dynamical fermions. The unique aspect of our work is the continuum extrapolation of all three 
renormalized quantities we have measured, chiral condensate, strange quark number susceptibility 
and Polyakov loop, which is especially difficult with Wilson fermions. The difficulty lies in the 
fact that the explicit chiral symmetry breaking of the Wilson formulation necessitates fine lattice 
spacings and consequently large lattice volumes. The continuum extrapolation has been carried out 
for all temperatures and by using four lattice spacings. 

The motivation for our work is twofold. First, it is desirable to obtain continuum QCD results 
from first principles which do not contain theoretically not fully justified operations such as the 
fourth root trick of rooted staggered fermions. The Wilson fcrmion formulation is theoretically 
sound and is known to be in the right universality class for QCD. Second, since a large body of 
simulation results for both zero and finite temperatures exist with the staggered formulation it is 
a worthwhile task to compare some of these continuum results with the continuum Wilson results. 
The expectation is of course that if the staggered formulation is also in the right universality class 
the continuum results will agree within errors. 

The Wilson formulation is more expensive than the staggered formulation hence lowering the 
light quark masses towards their physical value is much more difficult . In this work the pion mass 
has been set to 545 MeV, heavier than physical, but it is important to note that at arbitrary quark 
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Figure 3: Renormalizcd Polyakov loop. The top left panel shows the Wilson results at the four 
lattice spacings corresponding to four (3 values. The continuum extrapolated result is also shown 
by the solid band. The top right panel shows the staggered results in the fixed- N t approach also 
together with the continuum extrapolated result. The bottom panel compares the two continuum 
results. 
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masses the continuum results are universal. Hence provided the same renormalization prescription 
is used for two discretizations the results should still agree in the continuum for arbitrary quark 
masses. We have carried out staggered simulations using the same renormalization prescription as 
we did for the Wilson simulations and have compared the continuum extrapolated results. Nice 
agreement was found for all three quantities. 

In the future we plan to lower the light quark masses towards their physical value. With an 
unimproved fcrmion action this would be almost hopeless however the use of stout improvement 
makes the fluctuation of low lying Dirac eigenvalues much less allowing for a smaller pion mass at 
a given volume. 

We also plan to compare the continuum extrapolated staggered and Wilson results with a formu- 
lation which combines the key advantages of the two: chiral symmetry (staggered) and theoretical 
soundness (Wilson). Such a lattice chiral discretization is the overlap formalism. Results at two 
lattice spacings are already available in the fixed- N t approach [33] ■ Once completed with a third 
or fourth lattice spacing a continuum extrapolation will be possible allowing for the comparison of 
three continuum results from three discretizations of QCD. 
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